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The adiabatic approximation in time-dependent density functional theory (TDDFT) yields reliable 
excitation spectra with great efficiency in many cases, but fundamentally fails for states of double- 
excitation character. We discuss how double-excitations are at the root of some of the most challenging 
problems for TDDFT today. We then present new results for (i) the calculation of autoionizing reso- 
nances in the helium atom, (ii) understanding the nature of the double excitations appearing in the 
quadratic response function, and (iii) retiieving double-excitations through a real-time semiclassical 
approach to correlation in a model quantum dot. 
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I. INTRODUCTION 

There is no question that time-dependent density 
functional theory (TDDFT) has greatly impacted calcu- 
lations of excitations and spectra of a wide range of sys- 
tems, from atoms and molecules, to biomolecules and 
solids HH. Its successes have encouraged bold and ex- 
citing applications to study systems as complex as pho- 
tosynthetic processes in biomolecules, coupled electron- 
ion dynamics after photoexcitation, molecular trans- 
port, e.g. Refs. 10-13] • The usual approximations used 
for the exchange-correlation (xc) potential in TDDFT cal- 
culations however perform poorly in a number of sit- 
uations particularly relevant for some of these appli- 
cations, e.g. states of double-excitation character 0], 
long-range charge-transfer excitations IS-iSJ, conical in- 
tersections 0, llOll , polarizabilities of long-chain poly- 
mers lITTI] , optical response in solids fT2L[l3ll . Improved 
functionals, modeled from first-principles, have been 
developed, and are beginning to be used, to treat some 
of these. In this paper, we focus on the problem of 
double-excitations in TDDFT, for which recently much 
progress has been made. 

The term "double-excitation" is a short-hand for 
"state of double-excitation character". In a non- 
interacting picture, such as Kohn-Sham (KS), a double- 
excitation means one in which two electrons have been 
promoted out of orbitals occupied in the ground-state, 
to two virtual orbitals, forming a "doubly-excited" 
Slater determinant. The picture of placing electrons in 
single-particle orbitals however does not apply to true 
interacting states. Instead one may expand any inter- 
acting state as a linear combination of Slater determi- 
nants, say eigenstates of some one-body Hamiltonian, 
such as the KS Hamiltonian, and then a state of double- 
excitation character is one which has a significant pro- 
portion of a doubly-excited Slater determinant. Clearly, 
the exact value of this proportion depends on the level 
of theory used for the groimd-state reference e.g. it will 
be different in Hartree-Fock than in TDDFT. 

Given the important role of correlation in states of 



double-excitation character, the question of how these 
appear in a single-particle-based theory has both fim- 
damental and practical interest. We shall begin by re- 
viewing the status of linear-response TDDFT in this re- 
gard: one must go beyond the ubiquitous adiabatic ap- 
proximation to capture these. We conclude the introduc- 
tion by discussing systems where double-excitations are 
particularly important, and we shall see that some of 
these cases are due to the peculiarity of the KS single- 
Slater determinant in the ground-state (e.g. in certain 
long-range charge-transfer states). Then in Section |ll] 
we consider what happens when the double-excitation 
lies in the continuum, and test a recently developed ker- 
nel approximation to describe the resulting autoioniz- 
ing resonance in the He atom. In Section HUl we study 
whether adiabatic kernels can be redeemed for double- 
excitations by going to quadratic response theory. Sec- 
tion HV] turns to a new approach that was recently pro- 
posed for general many-electron dynamics, that uses 
semiclassical dynamics to evaluate electron correlation. 
Here we test it on a model system to see whether it cap- 
tures double-excitations. 



A. Double excitations in linear-response TDDFT 

Excitation spectra can be obtained in two ways from 
TDDFT. In one, a weak perturbation is applied to the 
KS system in its ground-state, and the dynamics in real- 
time of, e.g. the dipole moment is Fourier transformed 
to reveal peaks at the excitation frequencies of the sys- 
tem, whose strengths indicate the oscillator strength. 
More often, a formulation directly in the frequency- 
domain is used, in which two steps are performed: 
First, the KS orbital energy differences between occu- 
pied (i) and unoccupied (a) orbitals, ujs = — are 
computed. Second, these frequencies are corrected to- 
wards the true excitations thr oug h solution of a general- 
ized eigenvalue problem fT3,ll5|] , utilizing the Hartree- 
exchange-correlation kernel, /Hxc["-o](r, r', cj) = l/|r — 
r'l + f xcinolir, r' ,Ld), a functional of the ground-state 
density no(r). 



Fundamentally, the origin of the linear response for- 
malism is the Dyson-like equation that links the density- 
density response function of the non-interacting KS sys- 
tem, Xs, with that of the true system, x- 



= Xs(w) +Xs(w)*/hxc(w)*x(w) 



(1) 



where we use the shorthand * to indicate the integral, 
Xs('^)*/hxc(w) = / d^ri Xs(r,ri,w)/Hxc(ri,r',w) think- 
ing of Xi Xsj /hxc etc as infinite-dimensional matrices in 
r, r', each element of which is a function of lu. The inter- 
acting density-density response function x[no]{j^, ^' ,t — 
t') = 5n{r,t)/Sv{Y't')\n=no measures the response in the 
density n(r, t) to a perturbing external potential v(r, t). 
In the frequency domain. 



X(r,r',w) = ^ 



(^'o|n(r)|*/)(*/|n(r')|«'o> 



(vl/o|n(r')|vI//)(*,|n(r)|*o) 



^0^ 



(2) 



where / labels the interacting excited-states, and uji = 
Ei — Eq is their frequency relative to the ground-state. 
Similar expressions hold for the KS system, substitut- 
ing the interacting wavefunctions above with KS single- 
Slater-determinants. 

In almost all calculations, an adiabatic approxima- 
tion is made for the exchange-correlation (xc) kernel, 
i.e. one that is frequency-independent, corresponding 
to an xc potential that depends instantaneously on the 
density. But it is known that the exact kernel is non- 
local in time, reflecting the xc potential's dependence 
on the history of the density. It is perhaps surpris- 
ing that the adiabatic approximation works as well as 
it does, given that even weak excitations of a system 
lead it out of a ground-state (even if its density is that 
of a groimd-state). One of the reasons for its success 
for general excitation spectra is that the KS excitations 
themselves are often themselves good zeroth-order ap- 
proximations. But this reason cannot apply to double- 
excitations, since double-excitations are absent in lin- 
ear response of the KS system: to excite two electrons 
of a non-interacting system two photons would be re- 
quired, beyond linear response. Only single-excitations 
of the KS system are available for an adiabatic kernel 
to mix. Indeed, if we consider the KS linear density re- 
sponse function, the numerator of Eq. |2]with KS wave- 
functions contains ($o|'^(r)|$/), which vanishes if the 
excited determinant $/ differs from the ground-state $o 
by more than one orbital. The one-body operator h{r) 
cannot connect states that differ by more than one or- 
bital. The true response fimction, on the other hand, re- 
tains poles at the true excitations which are mixtures of 
single, double, and higher-electron-number excitations, 
as the numerator (^'o|n.(r)|4'7) remains finite due to the 
mixed nature of both and ^P/. Within the adiabatic 
approximation, x therefore contains more poles than Xs . 

Ref. iT6ll pointed out the need to go beyond the 
adiabatic approximation in capturing states of double- 



excitation character. Ref. derived a frequency- 
dependent kernel, motivated by first-principles, to be 
applied within the subspace of single KS excitations that 
mix strongly with the double-excitation of interest. For 
the case of one single-excitation, g = i — > a, coupled to 
one double-excitation: 



2tel/xcHI'?]=2[<z|/4K)k] 



\H. 



qD\ 



DD 



(3) 



to be applied within a dressed single-pole approxima- 
tion ("DSPA"), 



UJ =UJq+ 2[9|/xc(w)|q] 



(4) 



The Hamiltonian matrix elements in the dynamical 
correction (second term of Eq. are those of the 
true interacting Hamiltonian, taken between the sin- 
gle (q) and double (D) KS Slater determinants of inter- 
est, as indicated, and Hqq is the expectation value of 
the true Hamiltonian in the KS ground-state. This can 
be generalized to cases where several single-excitations 
and double-excitations strongly mix, within a "dressed 
Tamm-Dancoff " scheme (see, e.g. (TtLIT^ ). The kernel 
is to be applied as an a posteriori correction to a usual 
adiabatic calculation: first, one scans over the KS orbital 
energies to see if the sum of two of their frequencies lies 
near a single excitation frequency, and then applies this 
kernel just to that pair. 

Essentially the same formula results from derivations 
with different starting points: in Ref. [|19|] , it emerges as a 
polarization propagator correction to adiabatic TDDFT 
in a superoperator formalism, made more rigorous in 
Ref. HJ. Ref. iH utilized the Bethe-Salpeter equa- 
tion with a dynamically screened Coulomb interaction, 
while Ref. ||22f| extended the original approach of Ref. 0] 
by taking account of the coupling of the single-double 
pair with the entire KS spectrum via the common en- 
ergy denominator approximation (CEDA). 



B. When are double excitations important? 

Even the low-lying spectra of some molecules are 
interspersed with states of double-excitation character, 
but we will argue that they also lie at the root of sev- 
eral significant challenges approximate TDDFT faces for 
spectra and photo-dynamics. Although not tradition- 
ally seen as a double-excitation problem, we will see 
that double-excitations haunt the difficulty in describ- 
ing conical intersections and certain long-range charge- 
transfer states. 

Molecular spectra First, double-excitations in their own 
right are prominent in the low-lying spectra of many 
conjugated polymers. A famous case is the class of 
polyenes (see Ref. llTll for many references). For ex- 
ample, in butadiene the HOMO to (LUMO+1) and 
(HOMO-1) to LUMO excitations are near-degenerate 
with a double-excitation of the HOMO to LUMO. If one 
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runs an adiabatic calculation and simply assigns the en- 
ergies according to an expected ordering, one obtains 
7.02eV for the vertical excitation from a B3LYP calcu- 
lation (similar with other hybrid functionals), in a 6- 
311G(d,p) basis set, while a CASPT2 calculation yields 
6.27eV. By using different basis sets, a more accurate 
value can appear, but rather fortuitously, since the state 
obtained in adiabatic B3LYP has more of a Rydberg, char- 
acter, rather than double character ll23ll . In Ref. fl^ , the 
dressing Eq. Q, generalized to a subspace of two KS 
single excitations instead of one, was applied, yielding 
6.28eV. Similar successes were computed for hexatriene, 
and also for 0-0 excitations. This system was later stud- 
ied in detail in Ref. |T^. analyzing more fully aspects 
such as self-consistent treatment of the kernel, and use 
of KS versus Hartree-Fock orbitals in the dressing. Fur- 
ther, in Ref. i23l , excited state geometries were success- 
fully computed in this way. Most recently, an extensive 
study of Eq. (|3) was performed in Ref. [25] to low-lying 
states of 28 organic molecules. 

Charge-transfer excitations It is well known that long- 
range charge-transfer excitations are severely underes- 
timated with the usual approximations of TIDDFT. The 
usual argument to explain this is that the TDDFT cor- 
rection to the bare KS orbital energy difference van- 
ishes because the occupied and unoccupied orbitals, one 
being located on the donor and the other on the ac- 
ceptor, have negligible overlap as the donor-acceptor 
distance increases d H m. The TDDFT prediction 
then reduces to the bare KS orbital energy difference, 
cl (acceptor) — e^f (donor) where the subscripts H and L 
refer to HOMO and LUMO, respectively. This typically 
leads to an underestimation, because in usual approxi- 
mations |e^f I underestimates the true ionization energy, 
while the lowest unoccupied molecular orbital (LUMO), 
el/ lacks relaxation contributions to the electron affinity. 
The last few years have seen many methods to correct 
the underestimation of CT excitations, e.g. Refs. ||27| - 
I29II : most modify the ground-state functional to correct 
the approximate KS HOMO's underestimation of /, and 
mix in some degree of Hartree-Fock, and most, but not 
all I2^l29ll determine this mixing via at least one empir- 
ical parameter. 

But the argument above only applies to the case 
where the donor and acceptor are closed-shell species; 
instead, if we are interested in charge-transfer between 
open-shell species (e.g. in something like LiH), the 
HOMO and LUMO are delocalized over the whole 
molecule. This is the case for the exact ground-state KS 
potential, as well as for semi-local approximations jSOll . 
The exact KS potential has a peak and a step in the bond- 
ing region, that has exactly the size to realign the atomic 
HOMO'S of the two fragments [30-32] (see also Fig.^. 
As a result the molecular HOMO and LUMO are delo- 
calized over the whole molecule. The HOMO-LUMO 
energy difference goes as the tunnel splitting between 
the two fragments, vanishing as the molecule is pulled 
apart; therefore every excitation out of the KS HOMO 
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FIG. 1: Excitations out of a heteroatomic molecule composed 
of open-shell fragments at long-range (eg. "stretched" LiH). 
Blue denotes orbitals occupied in the ground-state; in this 
model, we show one electron on each "atom" (the inner elec- 
trons play only a secondary role). The top panel shows a 
model of the possible excitations of the true system: on the 
left is a local excitation on one atom, on the right are shown 
two lowest charge-transfer excitations. The bottom left panel 
shows the corresponding KS potential, displaying a step and 
a peak, as discussed in the text. The groimd-state KS is the 
doubly-occupied bonding orbital; any single excitation out of 
here is near-degenerate with a double-excitation where the 
other electron occupying the bonding orbital transits to the 
near-degenerate antibonding LUMO. 



is near-degenerate with a double-excitation where a 
second electron goes from the HOMO to the LUMO 
(at almost zero KS cost). This KS double-excitation 
is crucial to capture the correct nature of the true ex- 
citations as otherwise we are left with half-electrons 
on each fragment, physically impossible in the dis- 
sociated limit. This strongly affects the kernel struc- 
ture, imposing a severe frequency-dependence for all 
excitations, charge-transfer and local, for heteroatomic 
molecules cornpqsed of open-shell fragments at large 
separation ||33|,IM1 . 

Coupled electron-nuclear dynamics The importance of 
double-excitations for coupled electron-nuclear dynam- 
ics was highlighted in Ref. [10]: even when the verti- 
cal excitation does not contain much double-excitation 
character, the propensity for curve-crossing requires ac- 
curate double-excitation description for accurate global 
potential energy surfaces. The same paper pointed out 
the difficulties TDDFT has with obtaining conical inter- 
sections: in one example there, the TDDFT dramatically 
exaggerated the shape of the conical intersection, while 
in another, its dimensionality was wrong, producing a 
seam rather than a point. Although a primary task is to 
correct the ground-state surface, the problem of double- 
excitations is likely very relevant arotmd the conical in- 
tersections due to the near-degeneracy. 

Autoionizing resonances In the next section we dis- 
cuss the case when the excitation energy of a double- 
excitation is larger than the ionization energy of the sys- 
tem. In this case, an autoionizing resonance results. 
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II. AUTOIONIZING RESONANCES IN THE HE ATOM 

Photoionization has a rich history in quantum me- 
chanics, with the photoelectric effect playing a pivotal 
role in establishing the dual nature of light, and contin- 
ues to be a valuable tool in analyzing atoms, molecules, 
and solids. The photospectra above the ionization 
threshold are characterized by autoionizing resonances, 
whereby bound state excitations interact with those into 
the continuum. Many theoretical methodsfH-SSl have 
been developed in order to predict these resonances. 

As TDDFT is both accurate and relatively computa- 
tionally inexpensive, it stands as a useful candidate for 
studying these excitations. For autoionizing resonances 
where a single excitation frequency (e.g. core Ryd- 
berg) lies in the continuum, TDDFT has been shown to 
work well|41-47]. However it was noted ll4TI,l43ll that res- 
onances arising from bound double excitations are com- 
pletely missing in the adiabatic approximation (as we 
might expect from section ITB)| . 

Ref . [43 derived a frequency-dependent kernel which 
allows TDDFT to predict bound-double autoionizing 
resonances. Below we review this derivation before test- 
ing it on the Helium atom. 



Frequency-dependent kernel for autoionizing 
double-excitations 



A. 



Fano's pioneering work on photoionization 14911 con- 
sidered a zeroth-order unperturbed system with a 
boimd state degenerate with that of a continuum 
state ^E- The unperturbed system differs from the true 
system by the coupling term Vcpi, and we define the ma- 
trix element Ve = {^E\Vcpi\^b) between the two states. 
The transition probability, for some transition opera- 
tor T (e.g. the dipole), between an initial state \i) and 
the mixed state with energy in the continuum was then 
found to be 



where 



E- Er 



r/2 



and the energy of the resonance is shifted to 



Er=Eb+V / dE- 



E-E' 



The asymmetry of the lineshape is given by 

^ ($b|r|z) + P / cLE'Ve' {<^e' \f\i)/{E - E') 
nVE{-l>E\f\i) 



(5) 



(6) 



(7) 



(8) 



while the width of the resonance is given by 

r = 27r|($B|14pi|*b)l'-2^l^ij|' (9) 

In Ref. this formalism is applied to the KS sys- 
tem, where $6 and $ e are now interpreted as KS wave- 
functions, and the full Hamiltonian is related to the KS 
Hamiltonian via 



Vcpl = Vcc ~ Vh- Vxc (10) 

By comparing Eq. ^ with the TDDFT linear response 
equations, a frequency-dependent kernel, valid in the 
region near the resonance, was derived 



/hxc(w) = Xs ^ - Xs + 



r{r/2 + i{uj - ujr)) 



(11) 

In Eq. [HI g = 1 is required by the facts that both the 
KS system and the true system respect the Thomas- 
Reiche-Kuhn sum rule and that the double-excitation 
does not contribute to the oscillator strength of the KS 
system [48]. This kernel can be used to 'dress' the ab- 
sorption spectra found via the adiabatic approximation 
(A A), a^{uj), giving an absorption spectra 



r/2)2 



(12) 



which contains both the AA excitations and the autoion- 
ization resonances. 

A few points are worth noting. Fano's zeroth order 
picture is assumed to take care of all interactions except 
for the resonance one. When we approximate this by 
the non-interacting KS one, we can therefore only ex- 
pect the results to be accurate in the limit of weak in- 
teraction, where the dominant interaction is the reso- 
nant coupling of the boimd double-excitation with the 
continuum states. When the kernel is included on top 
of an adiabatic one, that includes also mixing of non- 
resonance single excitations, it is done in an a posteriori 
way, i.e. non-self-consistently, and expressions for the 
width are unaltered, do not include any adiabatic cor- 
rection. It is also worth noting that the derivation con- 
siders an isolated resonance: just one discrete state and 
one continuum. 



B. Application: Is^ 2s^ resonance in the He atom 
In order to test the accuracy of this prescription, we 



studied the l; 



2s autoionization resonance of He- 



lium. This is the lowest double-excitation in the He 
atom, and it lies in the continuum: Experimentally this 
excitation occurs at a frequency of w = 57.82 eV while 
the ionization threshold is / = 24.59 eV. In this case the 
wavefunctions are given by 

$b(r,r') - <^2s(r)</'2s(r') (13) 
$ij(r,r') = -i=(0is(r)0B(r') + </>is(r')0£(r))(14) 
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where 0i5(r) and 02s(r) are bound KS orbitals while 
(/)B(r) is the energy-normahzed continuum state with 
energy E. For various xc functionals, the bound state 
orbitals were calculated using OCTOPUS[50] while the 
unbound state was found using an RK4 integrator on 
the KS potential. 

In order to produce a continuum state with the cor- 
rect asymptotic form, the KS potential should decay 
like — 1/r. Two fimctionals which meet this condi- 
tion for Helium are the local density approximation 
with the Perdew-Zungerfsill self-interaction correction 
(LDA-SIC) and the exact exchange functional within the 
optimized effective potential method (OEP EXX). 

In table IB we show the results for the autoionization 
width found using Eq. (|9j. As can be seen the widths 
for the transition to the 2s^ state are too low Closer in- 
spection of the bound-state electronic structure reveals 
that there may be significant mixing between the 2s^ 
state and the 2p^ state. So we diagonalize the full 
Hamiltonian in this two-by-two subspace in order to in- 
clude the effect of this configuration interaction mixing. 
From this diagonalization, we take the state dominated 
by 2s^ wavefunction (about 70% in our cases). The au- 
toionization width for this mixed state (denoted "two- 
configuration") is also shown in TableHwhere it indeed 
improves the pure 2s^ state results, but is still roughly 
a factor of 40% too small. The LDA-SIC functional per- 
formed better than exact-exchange, this may be due to 
correlation being strongest in the core region which con- 
tributes most to the integral of Eq. ||9]l. The width for the 
2P^ excitation is too small to be measured experimen- 
tally, however we can compare to other theoretical cal- 
culations. After diagonalization, the width does become 
much smaller compared to the pure state, as it should, 
however the value is an order of magnitude smaller than 
that of Ref.lH 



TABLE I: Autoionization widths for Helium as calculated by 
various functionals. 



In conclusion, we tested the formalism of Ref. 148 
to include double-excitation autoionization resonances 
within TDDFT for the Helium atom. The results, while 
not outlandish, are disappointing when compared to 
simpler wavefunction methodsfsS] for this resonance in 
Helium. This is probably due to the fact that the KS 
system in this case does not make a good zeroth-order 
picture on which to build a Fano formalism: i.e. the 
assumption of weak interaction, aside from the reso- 
nance coupling, discussed in Sec. Ill Al do not hold. It 
would be interesting to compute the shift in the res- 
onance position (Eq. lO) within the TDDFT prescrip- 
tion; this is expected to be large, e.g. comparing the 
KS double (w = 40.6 eV in exact TDDFTIS^) to the ex- 
perimental resonance (cj = 57.82 eV). Two approximate 
functionals were tested, LDA-SIC and exact exchange, 
with the results suggesting that correlation improves 
the description. Functionals missing the — 1/r tail in 
their KS potentials can provide continuum states accu- 
rate within a core regionfs^ IstII (although missing the 
correct asymptotic behavior) and so may still be used 
with this method. However they are unlikely to sig- 



state 



Method 



r(ev) 



2s' 



OEP EXX pure 

OEP EXX two-configuration 

LDA-SIC pure 

LDA-SIC two-configuration 

MCHF+b-spline" 

Experiment'' 



0.0604 
0.0741 
0.0670 
0.0836 
0.1529 

0.138 ±0.015 



2p' 



OEP EXX pure 0.0139 

OEP EXX two-configuration 0.0003 

LDA-SIC pure 0.0170 

LDA-SIC two-configuration 0.0004 

MCHF+b-spline" 0.0055 

Experiment - 

"From Ref. 152] ''From Ref. Q 



nificantly improve the results for the reasons discussed 
above. 

Hellgren and von Barth derived the Fano lineshape 
formula in terms of linear response quantities, that is 
exact within the adiabatic approximation ji^ . It can- 
not apply to the case of a bound double-excitation. On 
the other hand, our expression does apply but it is a 
lower-order approximation, only expected to be accu- 
rate in the limit of a narrow isolated resonance in the 
weak interaction limit. For molecules and larger sys- 
tems, TDDFT is often the only available technique to 
calculate excitations and this formalism is the only avail- 
able method to treat doubly-excited autoionization res- 
onances within TDDFT. For such larger systems, we ex- 
pect this approach might provide useful results. 



III. ADIABATIC QUADRATIC RESPONSE: DOUBLE 
VISION? 



Given that to excite two electrons in a non-interacting 
systems, two photons are required, one may ask 
whether adiabatic TDDFT in nonlinear response, espe- 
cially quadratic response, yields double-excitations. In 
this section, we work directly with the second-order KS 
and interacting response functions to investigate this 
question. We will first consider the form of the non- 
interacting quadratic response function and briefly re- 
view TDDFT response theory, before turning specifi- 
cally to the question of the double-excitations. 

Applying a perturbation to a system Sv{r,t) initially 
in its groimd state tiq (r), we may expand the response of 
the density in orders of Sv{r, t), where n(r, t) ~ no(i") + 
ni(r, t) + n2(r, + .... We define corresponding response 



5 



functions: 



X(r,t,r',t') = 



6n[v] (r, t) 



Sv{r',t') 

eit-t'){'^o\[nHir,t),hHir',t'mo) 



(15) 

as the linear-response density-density response function 
whose frequency-domain version appeared earlier in 
Eqs. lUIll, and the second-order response function 



x(2'(r,i,r',t',r",i") 



1 



S^n[v]{r, t) 



5v{v' ,t')5v{Y" ,t") 



V=Vo 



-{{^o\[[nH (r, t) ,nHir',t% hn (r" , )] | *o) 

+pcrmutations of (r', t') ^ (r", t")} e{t - t')9{t - t") 

(16) 

where n^f is the density operator in the Heisenberg 
picture and the expressions in terms of the commu- 
tators follow from time-dependent perturbation the- 
ory (H HI. 

In the following we often will drop explicit spatial- 
dependences, and think of density as a vector. Kernels 
and response functions would then be represented as 
matrices and tensors with the symbol * signaling con- 
traction. 



Slater determinant. Expanding out the commutator in 
Eq. | [T6l l and inserting the identity in the form of com- 
pleteness relations of the Slater-determinant basis, we 
obtain terms of the form 

^5](0s|n(r)|/s)(/s|n(r')|Js)(Js|n(r")|0s) (17) 
Is Js 

where the subscript S indicates it's a KS state. Exam- 
ining the first bracket, we see that Is must be either the 
ground-state or a single-excitation, since the one-body 
density operator can only connect determinants differ- 
ing by at most one orbital. Likewise, for a non-zero 
third bracket, can only be the ground-state or a single- 
excitation. Therefore, no double-excitations contribute 
to the second-order response function. At second-order, 
double-excitations are reached in the dynamics of a non- 
interacting system but cannot contribute to the one- 
body density response. By similar arguments, the low- 
est order density-response that double-excitations can 
appear is at third-order. That the KS second-order re- 
sponse fimction does not contain double-excitations be- 
gins to dash any hope that adiabatic TDDFT will yield 
accurate excitations. It will turn out that we do "see dou- 
ble", but the vision is blurry: they are simply sums of 
linear-response corrected single-excitations, quite blind 
to any nearby single excitation. Before turning to a 
closer investigation of this, we briefly review the struc- 
ture of TDDFT response theory. 



A. The non-interacting quadratic response function 

First, we ask, what is happening at the level of the 
non-interacting KS response functions? In the above 
equations, the functional derivatives with respect to v 
are replaced by those with respect to the corresponding 
Vs and the ground-state in Eqs. l|15H16t is then a single 



B. Non-linear response in TDDFT 

The fact that the time-dependent KS system repro- 
duces the true system's density response to the cor- 
responding KS perturbing potential (5ws(r, i]/ leads to 
Eq. ^ for the linear response function, and lISSll 



dsdf-y I ds' d^y'xs\r,t,y,s,y',s 



Sv{r',t') 



no Hr",t") 



+ dsd'yxsivAy.s) ds'd'y' ds" d^y VcSKy, s, y', s', y", s")x(y', s', r', i')x(y", s", r", i") 



+ J dsd^yxsir,t,y,s) J ds' d'y' U^c.[no]{y, s,y' , s')x^^Hy' , s' ,r' ,t' ,r" ,t") 



(18) 



S-'Excln] 



for the quadratic response function. Here Xs^^ denotes is the dynamical second-order xc kernel. For simplic- 

the second-order KS response function, discussed in ity we assume the state is spin-saturated. In the adi- 
SecHlIAland 

5xcM(r,i,r',t',r",t") 



d'^Vxc[n]ir,t) 



Sn{r',t')Sn{r",t") 



abatic approximation, g^^{r,r',r") = sn{r)Sn(r')5n{r") 
with i?xc[«] a ground-state energy functional. Making 
a Fourier transform with respect to t — t' and t — t" we 



(19) 
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KS determinants and KS energies, one finds that, for uj 
not equal to a KS transition frequency or its negative. 



XsH-E^^ (25) 



where q labels a transition from occupied orbital i to an 
unoccupied orbital a and the matrix 



As,,(r,r') =4w,<^,;(r)<^a(r)</)a(r')04r') (26) 



where Wg = Ca — and we take the KS orbitals to be real. 
We note that there is a factor of 2 in As,q(r, r') due to the 
assumed spin-saturation. Also note that if only forward 
transitions are kept, this amounts to what is often called 
the Tamm-Dancoff approximation: 



aTD 

Xs(c.)=E^^ (27) 
q ^ 



with^j;^ = 2Uv)^a{v)K{r^')Uv'). 

Returning to Eq. |24l for an adiabatic approximation, 

hxci^) = /hxc we may write 



/ J p A \ 



obtain the first and second-order density responses as 

ni{u!) = x(aj) * (5w(a;) 

= Xsi^j) Sv{oj) + Xs('^)*/Hxc(w)*ni(a;) 

(20) 

and 

712(0;) ~ ^'^'x^'^\'^7'^ '-^') * ^ '^') 

+ ^Xs(w) * gti^cioJ,uj - uj')-k J duj'ni(uj')ni{u} - uj') 

+Xs(a;) * /hxc(w) *n2(cj) (21) 
respectively. 

Ref . Is^ pointed out a very interesting structure that 
the TDDFT response equations have. At any order i, 

ni{uj) = Mi{ui) + Xs('^) -* /hxc('^) *ni(w) (22) 

where A/, depends on lower-order density-response (and 
response-functions up to ith order). The last term on the 
right of Eq. |22]has the same structure for all orders. If 
we define the operator 

L(w) = 1-Xs(w)*/hxc (23) 

then 

L{uj)i.n,{uj)=M,{uj) (24) 

In the next section we work out the effect of this oper- 
ator in the adiabatic approximation, by studying linear 
response. This will be useful to us when we finally ask 
about doubles in the quadratic response function. 



Xs(c.)^/h1cc-E 2 (28) 

q 1 

whereAs*/4c = 4a;g0i(r)(/)a(r) / /4c(ri, r')'/'»(ri)(/)a(ri)fi3ri. 



This means the effect of the operator on the left-hand- 
C. Adiabatic approximation ^ide of Eq. El is to zero out poles of the fimction it is 

operating on that lie at KS single excitations, and re- 
Let us find the inverse of the operator appearing on place them with linear-response-corrected ones, within 
the left-hand-side of Eq.|24l First, evaluating Eq.|2]with the adiabatic approximation. In particular. 
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(where Ilg indicates a product over q) It is instructive andl20l 
to first consider linear response (z = 1 in Eq. l24t , and _^ 
zoom in on only one excitation, with coupling to all oth- ni{uj) = L (uj) * Xsi^^) * Sv{u}) 
ers considered insignificant. Then, from Eqs.|28l EH |29l ((w^ - w^)! - (As ★ /^c)?) "^ * ^s,? * '5«(^) 

(30) 
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Poles of ni (w) are thus indicated by zeroes of the first 
term in brackets: when (oj^ — a;^)(5(r — r') — (As*/^xc)g = 
0. Integrating this over r2 and reahzing the second term 
gives zero unless r' = r, we find we have effectively 
derived the "small-matrix approximation" ll60l,l6lll 

(31) 

Had we kept only forward terms, and proceeded in a 
similar manner, the single-pole approximation (Eq.|4]but 
with a frequency-independent right-hand-side) would 
have resulted. The essential point is that the effect of 
L^^ is to shift the KS pole towards the linear-response 



corrected single excitation. 

Now this operator appears at all orders of TDDFT 
response, yielding critical corrections to the non- 
interacting response functions at all orders. We now are 
ready to investigate the question of double-excitations 
in quadratic response in the adiabatic approx. 



D. The quadratic density-response in the adiabatic 
approximation 

We have from Eq. l(2T)l and definition l|23] | that 



n2{uj) = - (^L "'"(w)* J dw'xs^'(w — w'jw) *(5vs(w — w') *(5us(i^') + ^i^)* Xs('^)*gHxc* J duj'ni{uj')ni{uj — uj') 



We wish to investigate specifically, the question of how 
double-excitations would appear in the poles of n2(w). 
From Sec. 1111 CI we know the effect of the operator 
L^^{ll!) acting on a function with a pole at a KS single 



excitation, is to shift that pole to its linear-response cor- 
rected value, in the adiabatic approximation, (Eqs. | |29l 
and l|30t . Let us denote these values as flj. We shall now 
study in detail the second term in Eq. 1111 Dl which is 



So until the last integral over lj', only poles at sin- response function can be written as 
gle KS excitations corrected by linear-response adiabatic 
TDDFT appear. We now turn to this last integral to see 

what it gives. First, we write it as x^{w) ~ ( — — ) (35) 

^ \uj — Qj + iO+ UJ + Qi + iO+ / 

d^yd^y' J duj'xir, y, uj')x{r' , y', uj-uj')G{y, y' ,uj' ,uj-uj') 

(34) which contains the same number of poles as Xsi^^)/ in 

where G'(y, y', w', w — a;') = wi(y,a;')ui(y',aj — a;'). Af- contrast to Eq. Inputting this adiabatic linear- 

ter applying an adiabatic kernel in Eq. the linear response fimction into Eq. (|34] |, we find: 



/ duj^yyAjAj ( — i i — 

Doing the integrals, we obtain 



1 1 



cj' — cj — rij uj' — (jj + Qj 



G{u',u-oj') (36) 



^^'l^l^Aj Aj y _____ + _____ 



G{-ni,Lo + ni) - Gjoj + nj, -Qj) Gjcj - Oj, nj) - G{~ni, Lo + ni) \ 
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This displays poles at sums and differences of the linear- 
response-corrected single-excitations. In particular, it 
contains poles at double-excitations, ±(17/ + Oj). These 
poles remain after being multiplied by the preceding 
terms in Eq. IIII Dl (which, as explained before, may can- 
cel/shift a pole at a single KS excitation). Notice that 
if we made a Tamm-Dancoff-like approximation for x 
with respect to cj' then Eq.l34l would be 



ij 

yielding 



w' — w — 17 J 



G(w',w-w') 
(38) 



^ ^^^^ G{ni,u~^i)-G{uj^nj,-nj) ^^^^ 



IJ 



uj — [VLj — VLj) 



That is, the double excitations vanish in this Tamm- 
Dancoff-like approximation. 

Our findings here are not entirely new: in Ref. m 
very similar conclusions were reached, but by quite a 
different method. To our knowledge, an analysis based 
directly on the response functions, as above, has not 
appeared in the literature before. The formalism in 
Ref. 16211 was based on the dynamics of a system of 
weakly coupled classical harmonic oscillators, in a co- 
ordinate system defined by the transition densities. The 
linear and non-linear response of this system was shown 
to correspond to that of a real electronic system via 
adiabatic TDDFT. Closed expressions for the first, sec- 
ond, and third optical polarizabilities were found; in 
particular the second-order polarizability was shown 
to contain poles at sums of linear-response corrected 



single-excitations, and that if a Tamm-Dancoff approx- 
imation was made, these poles vanish. We have shown 
very similar results staying within the usual language 
of TDDFT and response theory, directly from analyzing 
the response functions. The relation between the Tamm- 
Dancoff-like approximation we make here and that re- 
ferred to in f62| remains to be examined in more detail. 

So we have shown that the quadratic response in 
TDDFT within the adiabatic approximation is unable to 
describe double excitations. While the second-order re- 
sponse in adiabatic TDDFT does contains poles at the 
sum of linear-response corrected KS frequencies, it com- 
pletely misses the mixing between single and double KS 
states needed for an accurate description of double ex- 
citations. Moreover, using the Tamm-Dancoff approxi- 
mation in linear response makes even these poles disap- 
pear. These conclusions were reached by looking at the 
pole structure of the response equations directly. Using 
this approach, we were able to see that the second-order 
Kohn-Sham response function does not contain poles at 
the sum of the KS frequencies, which preventing TDDFT 
from introducing the mixing mentioned above. 

IV. DOUBLE-EXCITATIONS VIA SEMICLASSICAL 
DYNAMICS OF THE DENSITY-MATRIX 

Recently, time-dependent density-matrix functional 
theory (TDDMFT) has been explored as a possible rem- 
edy for many of the challenges of TDDFT [63, 64]. The 
idea is that including more information in the basic vari- 
able would likely somewhat relieve the job of xc fimc- 
tionals and lead to simpler fimctional approximations 
working better. The (spin-summed) one-body density- 
matrix is 



p(r',r,t) 



j dx2--dxN'i'*{r'(Ti,X2--XN,t)'i'{rai,X2--XN,t) 



(40) 



r 



(using X = (r,cr) as spatial-spin index with J dx ^ 
12 a I d^''')/ while the density is only the diagonal ele- 
ment, n{r,t) = p{r,r,t). For example, immediately we 
see that the kinetic energy is exactly given by the one- 
body density-matrix as T = — J d^r^V^ p{r'r,t)\r'=r, 
while only approximately known as a functional of the 
density alone; and, only the non-interacting part of the 
kinetic energy is directly calculated from the KS orbitals 
(as — / dx(j)* (a;)V^(/),;(a;)/2), while it is unknown how 
to extract the exact interacting kinetic energy from the 
KS system. TDDMFT works directly with the density- 
matrix of the interacting system so is not restricted to 
the single-Slater-determinant feature of the TDKS sys- 
tem. For this reason, (TD)DMFT would be especially 
attractive for strongly-correlated systems, for example 



dissociating molecules and it was recently shown that 
adiabatic approximations in TDDMFT are able to cap- 
ture bond-breaking and charge-transfer excitations in 
such systems ll63l . However, it was also shown that adi- 
abatic approximations within TDDMFT still cannot cap- 
ture double-excitations 



Very recently a semiclassical approach to correlation 
in TDDMFT has been proposed 116511 , which was argued 
to overcome several failures that adiabatic approxima- 
tions in either TDDFT or TDDMFT have. The appli- 
cations in mind involved real-time dynamics in non- 
perturbative fields, for example, electronic quantum 
control via attosecond lasers, or ionization processes. 
In these applications, memory-dependence, inclu ding 
initial-state dependence, is typically important ll66l |67|] , 
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but lacking in any adiabatic approximation. More 
severely, the issue of time-evolving occupation numbers 
becomes starkrngly relevant: typically, even when a sys- 
tem begins in a state which is well-approximated by a 
single-Slater determinant, it will evolve to one which 
fundamentally involves more than one Slater determi- 
nant (eg. in quantum control of He from Is^ — ^ ls2p, 
or in ionization) Is^ . Although impossible when evolv- 
ing with one-body Hamiltonians such as in TDKS (thus 
making the job of the exact xc potential very difficult, as 
well as observables to extract information from the KS 
system), in principle TDDMFT can change occupation 
numbers. However it was recently proven that adiabatic 
approximations in TDDMFT cannot I6a-[70n . The semi- 
classical correlation approach of Ref. | |65|] incorporates 
memory, including initial-state dependence, and does 
lead to time-evolving occupation numbers, as has been 
demonstrated on model systems (7l|] . We now ask, does 
it capture double-excitations accurately? We shall use a 
model two-electron system to investigate this question, 
but first will review the method. Unlike the previous 

I 



It is convenient to decompose this as 

P2{x[,x'2,xi,X2;t) = p{x'2,x2;t)p{x[,xi;t) 
- P{x'i,x2;t)p{x2,xi;t) 

+ P2cix[,x'2,Xi,X2]t) (44) 

where we identify the first term as the Hartree piece, the 
second as exchange, and the last as correlation. If p2c 
is set to zero, one obtains Hartree-Fock; it is this term 
that Ref. fbBl proposed to treat semiclassically in order 
to capture memory -dep endence and time-evolving oc- 
cupation numbers ll65ll . We shall shortly see that, in 
contrast to the adiabatic approximations of this term, 
its semiclassical treatment also approximately captures 
double-excitations . 

There are various different semiclassical formula- 
tions, and the one we will explore here is known as 
"Frozen Gaussian" propagation, proposed originally by 
Heller ll72ll . This can be expressed mathematically as 
a simplified version of the Heller-Herman-Kluk-Kay 
(HHKK) propagator ll73l [tSI where the wavefunction at 
time t: 

**(x) = |^|p^(x|q,p,)Cq,p,,e^^'/^qoPo|vI',) (45) 

where {qt,Pt} are classical phase-space trajectories 
in 6A^-dimensional phase-space, starting from initial 
points {qoiPo}/ and '^i is the initial state. In Eq. l45l 



sections in the paper, this operates in the real-time do- 
main, instead of the frequency-domain. 



A. Semiclassical correlation in density-matrix propagation 
The equation of motion of p is given by 

ip{r',r;t) = f — ^ + Woxt(r; t) + — - i'oxt(r') 1 p(r', r 
+ J dr2fcc{r' ,r,r2)p2{r' ,r2,r,r2;t) {' 

where /cc(r'r,r2) = l/|r-r2| - l/|r'-r2| and /92 is the 
second-order reduced density matrix defined by: 

/32(ri,r^,ri,r2;<) = ^ P2ir[ai,r'2cr2,riai,r2cr2;t) 

C"l ,0"2 

(42) 

and 



(43) 

I 

(x|qp) denotes the coherent state: 

1/4 

(x|qp) = Y[ (— ) e-^(^^-«^)'+*f^(^^-«^)/'' (46) 

where -jj is a chosen width parameter. St is the clas- 
sical action and C^ p t is a pre-factor based on the 
monodromy (stability) matrix. The pre-factor is time- 
consuming to compute, and scales cubically with the 
number of degrees of freedom, but in the Frozen Gaus- 
sian approximation, is set to unity. Ref. | |65| ] gave an 
expression for the second-order density-matrix within 
a Frozen-Gaussian approximation, that, furthermore, 
takes advantage of the fact that there is some phase- 
cancellation between and 5** in Eq. |43]so that the re- 
sulting phase-space integral is less oscillatory. 

The idea is to extract the correlation term from the 
semiclassical d5mamics and insert it as a driving term 
in Eq. ||4T] | [65]. That is, we compute the semiclassi- 
cal first-order and second-order reduced-density matri- 
ces from Frozen Gaussian d5mamics, placing them in 
Eq. l l44l l which is inverted to solve for pfc • the other 
terms of Eq. l|44|l and (|4T] | are given exactly in terms of 
the one-body density-matrix, we use only the semiclas- 
sical expression for p2c when driving Eq. l(4T|l . 

It is interesting however to first ask how well semi- 
classical dynamics alone does. That is, without cou- 
pling to the exactly-computed one-body Hartree and 



P2{x[,x'2,xi,X2;t) 



N{N - 1) 



dx^..dxN'^*{x[,x'^, x^..X]s; t)-^{xi,X2,r2a2,X:i..XN] t) 
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exchange-terms in Eq. ||4T] |, how would the semiclassical 
calculation alone predict the dynamics? In particular, do 
we obtain double-excitations, and if so, how well. 



B. Double excitations from semiclassical dynamics 

We consider the following one-dimensional model of 
a two-electron quantum dot: 



H 



(47) 



using a soft-Coulomb interaction between the elec- 
trons. Starting from an arbitrary initial state, the inter- 
acting dynamics in this Hamiltonian can be numerically 
solved exactly, and we will compare exact results with 
those computed by Frozen Gaussian dynamics. 

First, by considering a non-interacting reference, we 
can identify where single and double excitations lie, 
which, when interaction is turned on, will mix. The 
level sketch in the upper right of Fig. |2] shows this: in 
the ground-state, both electrons occupy the lowest level, 
and shown is a single excitation to the second-lowest 
excited orbital (left) and a double-excitation to the first 
excited orbital (right). The non-interacting energies of 
these two states are near-degenerate, and mix strongly 
to give roughly 50:50 single:double mixtures for the true 
interacting states. Due to quadratic symmetry of the 
Hamiltonian, these states do not appear in the dipole re- 
sponse of the system, hence we look at the quadrupole 
moment: 



q{t) = dx x'^n{x; t) 



(48) 



We start in an initial state quadratically 'kicked' from the 
grotmd state, that is 



Vl/,(xi,X2) = e^''(^'+^')*o(a;i,X2) 



(49) 



where -q is chosen large enough to sufficiently populate 
the states we are interested in without being too large 
leading to higher-order response effects. A value of 77 = 
0.01 was used in our calculation. The Frozen Gaussian 
integral of Eq. | [45l l is performed by Monte-Carlo inte- 
gration using trajectories based on importance sampling 
the initial state overlap with the coherent state, and a 
width of 7 = 1 was used. In the example discussed be- 
low, 120000 trajectories were used, with the quadrupole 
moment not changing significantly if more trajectories 
are used. These trajectories are then classically propa- 
gated forward in time using a standard leapfrog algo- 
rithm. A total time of T = 200au was performed with 
a timestep of At = 0.001, although the wavefunction is 
only constructed every 0.1 au as this is sufficient to see 
the frequencies we are interested in. 

In Fig. 121 we show the power spectrum of the Fourier 
transform of the quadrupole moment for a Frozen Gaus- 
sian calculation with the parameters given above. Also 




Peaks at mixed single and double 
excitations: — \ — 



Exact 


AEXX 


SC 


DSPA 


2.000 


1.87 


1.99 


2.000 


1.734 




1.58 


1.712 



FIG. 2: The Fourier transform of the quadrupole moment com- 
puted via semiclassical (Frozen Gaussian) dynamics in the 
Hooke's quantum dot. We focus on the region of the spectrum 
where the excitations shown in the level sketch to the upper 
right lie as explained in the text. The table shows the val- 
ues of the excitations computed exactly, using adiabatic exact- 
exchange (AEXX), semiclassical (SC), and the dressed correc- 
tion of Eq. d-d (DSPA). 



given in Fig. |2] is a table comparing the exact, adia- 
batic exact-exchange (AEXX), semiclassical (SC), and the 
dressed correction (DSPA) discussed in the introduction. 
As stated earlier, an adiabatic approximation cannot in- 
crease the number of poles in the KS response function: 
AEXX can only shift the KS single excitation, yielding a 
solitary peak in between the two exact f requencies lTSll . 
In contrast to this, the Frozen Gaussian semiclassical re- 
sults can be seen to give two peaks in the right region. 
Although one peak is lower than in the exact case, this 
error may be lessened when p2c is used to drive the 
density-matrix propagation, with the one-body terms 
treated exactly quantum-mechanically instead of semi- 
classically (Sec. IIV At . The cj = 2 peak arises from ex- 
citation in the center-of-mass coordinates where the re- 
duced Hamiltonian is harmonic, given that the Frozen 
gaussian method is exact for such systems, we would 
expect this peak to be very accurate. The DSPA works 
extremely well in this case, although in more compli- 
cated systems one must search the KS excitation spec- 
trum in search of nearby doubles, as explained in the 
introduction, whereas they will appear naturally in the 
spectrum in the semiclassical approach. 

This example demonstrates that semiclassical correla- 
tion does capture double-excitations approximately, un- 
like any adiabatic approach in TDDFT or TDDMFT. We 
stress that this is the result from semiclassical dynamics 
alone (within the frozen gaussian approximation); in fu- 
ture work we will investigate whether the coupling to 
the exact one-body terms of Eq. [41] improve the results 
for the excitations. 
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V. SUMMARY AND OUTLOOK 

TDDFT is in principle an exact theory based on a 
single-particle reference: the exact functionals extract 
from the non-interacting KS system, the exact excita- 
tions and dynamics of an interacting electronic sys- 
tem. As such, it is both fundamentally and practi- 
cally extremely interesting how these fimctionals must 
look when describing states of double-excitation char- 
acter, particularly in linear response where no double- 
excitations occur in a non-interacting reference. The 
work of Ref. |0] that shows the form of the exact xc 
kernel and models an approximate practical frequency- 
dependent kernel based on this, has recently drawn 
some interest, both from a theoretical and practical 
point of view. We have discussed here how double- 
excitations are at the root of some of the most difficult 
problems in TDDFT today: long-range charge-transfer 
excitations between open-shell fragments, and conical 
intersections. 

The paper then described three new results in 
three different approaches involving double-excitations. 
First, we have applied a recently proposed ap- 
proach 10] to autoionizing resonances arising from 
double-excitations, to compute the width of the 2s^ reso- 
nance in the He atom. Although the results are not very 
accurate, predicting a 40% too narrow resonance, this 
approach is the only available one for this kind of reso- 
nance in TDDFT today. For larger systems, where the al- 
ternative wavefunction-based methods are not feasible, 
the approach of Ref. jHI might still be useful, despite the 
weak-interaction assumption (that is likely responsible 
for the error in the He case). Second, we showed that use 
of the adiabatic approximation in quadratic response 
theory yields double-excitations which are the sums of 
linear-response-corrected single-excitations. Within a 



Tamm-Dancoff-like approximation, even these doubles 
disappear. We argued that the KS quadratic response 
fimction does not have poles at KS double-excitations, 
which is behind the reason that we do not see the truly 
mixed single and double excitations in the adiabatic 
TDDFT quadratic response fxmction. Although simi- 
lar conclusions were reached in Ref. (62l], our analysis 
here proceeds in a very different manner: here, we fol- 
low more traditional response theory within DFT, with- 
out introducing new formalism. Finally, we investi- 
gated whether and how accurately double-excitations 
appear in a recently proposed semiclassical approach 
to correlation. This approach was originally proposed 
for general real-time dynamics [6^ , based on propaga- 
tion of the one-body-density-matrix. The correlation- 
component of the second-order density matrix, that ap- 
pears in the equation of motion for the first-order one, is 
computed via semiclassical dynamics. Here we showed 
that running semiclassical dynamics on the whole sys- 
tem does approximately capture double-excitations. Fu- 
ture work includes retaining exact dynamics for the one- 
body terms, leaving semiclassics just for the correlation 
component, according to the original prescription, to see 
if the accuracy of the states of double-excitation charac- 
ter are improved. 

In conclusion, there are many fascinating things to be 
learnt and discovered about states of double-excitation 
character! We hope that the findings here, and in the ear- 
lier work reviewed in the introduction, will spur more 
interesting investigations, with both practical and fim- 
damental consequences. 

We acknowledge support from the National Science 
Foundation (CHE-0647913), the Cottrell Scholar Pro- 
gram of Research Corporation and the NIH-funded 
MARC program (SG and CC). 
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